// 版权所有2010 Go作者。保留所有权利。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

/*
	Bessel function of the first and second kinds of order one.
*/

// 原始C代码和下面的长注释是
// 来自FreeBSD的/usr/src/lib/msun/src/e_j1。c和
// 随此通知一同发出。go代码是原始C.
// 
// ==============================================================================================================================
// Sun Microsystems，Inc.版权所有。
// 
// 由Sun Microsystems，Inc.的企业SunPro开发。
// 使用、复制、修改和分发本
// 软件的权限是免费授予的，前提是保留本通知
// 的权利。
// /==================================================================================================================
// /
// /\uU ieee754_j1（x），\uU ieee754_y1（x）
// Bessel函数的第一类和第二类一阶。
// Method--j1（x）：
// 1。对于微小的x，我们使用j1（x）=x/2-x**3/16+x**5/384-。。。
// 2。将x减少到| x |，因为j1（x）=-j1（-x），对于x in（0,2）
// j1（x）=x/2+x*z*R0/S0，其中z=x*x；
// （精度：|j1/x-1/2-R0/S0 |<2**-61.51）
// x在（2，inf）
// j1（x）=sqrt（2/（pi*x））*（p1（x）*cos（x1 q1（x）*sin sin x1。最好计算sin（x1），cos（x1）cos（x1）cos（x1）=cos（x）cos（3pi/4）+sin（x）sin（3pi/4）
// /=1/sqrt（2）*（sin（x）-cos（x））
// （sin x1）=sin x cos（3pi/4）-cos（x）sin（3pi/4）
// sin（x）+-cos（x）=-cos（2x）/（sin（x）-+cos（x））
// 来计算更差的一个。）
// 
// 3个特殊情况
// j1（nan）=nan 
// j1（0）=0 
// j1（inf）=0 
// 
// 方法-y1（x）：
// 1。筛选出x<=0例：y1（0）=-inf，y1（x<0）=NaN 
// 2。对于x<2。自
// y1（x）=2/pi*（j1（x）*（ln（x/2）+Euler）-1/x-x/2+5/64*x**3-…）因此y1（x）-2/pi*j1（x）*ln（x）-1/x是一个奇数函数。我们使用以下函数来近似y1，
// y1（x）=x*U（z）/V（z）+（2/pi）*（j1（x）*ln（x）-1/x，z=x**2 
// x在[0,2]中（绝对误差小于2**-65.89）
// U（z）=U0[0]+U0[1]*z++U0[4]*z**4 
// V（z）=1+v0[0]*z+…+v0[4]*z**5 
// 注意：对于微小的x，1/x支配y1，因此
// y1（微小）=-2/pi/tiny，（选择微小的<2**-54）
// 3。对于x>=2。
// y1（x）=sqrt（2/（pi*x））*（p1（x）*sin（x1）+q1（x）*cos（x1））
// 其中x1=x-3*pi/4。最好用上述方法计算sin（x1），cos（x1）
// 。

// J1返回第一类的一阶贝塞尔函数。
// 
// 特殊情况为：
// J1（±Inf）=0 
// J1（NaN）=NaN 
func J1(x float64) float64 {
	const (
		TwoM27 = 1.0 / (1 << 27) // 2**-27 0x3E40000000000000000 
		Two129 = 1 << 129        // 2**129 0x480000000000000 
		S01 = 1.91537599538363460805e-02  // 0x3F939D0B12637E53 
		S02 = 1.85946785588630915560e-04  // 0x3F285F56B9CDF664 
		S03 = 1.17718464042623683263e-06  // 0x3EB3BFF8333F8498 
		S04 = 5.04636257076217042715e-09  // 0x3E35AC88C97DFF2C 
		S05 = 1.23542274426137913908e-11  // 0x3DAB2CFCFBFB97ED8 
	)
	// 特殊情况
	switch {
	case IsNaN(x):
		return x
	case IsInf(x, 0) || x == 0:
		return 0
	}

	sign := false
	if x < 0 {
		x = -x
		sign = true
	}
	if x >= 2 {
		s, c := Sincos(x)
		ss := -s - c
		cc := s - c

		// 
		if x < MaxFloat64/2 {
			z := Cos(x + x)
			if s*c > 0 {
				cc = z / ss
			} else {
				ss = z / cc
			}
		}


		var z float64
		if x > Two129 {
			z = (1 / SqrtPi) * cc / Sqrt(x)
		} else {
			u := pone(x)
			v := qone(x)
			z = (1 / SqrtPi) * (u*cc - v*ss) / Sqrt(x)
		}
		if sign {
			return -z
		}
		return z
	}
	if x < TwoM27 { // /|x |<2**-27 
		return 0.5 * x // =0必需的
	}
	z := x * x
	r := z * (R00 + z*(R01+z*(R02+z*R03)))
	s := 1.0 + z*(S01+z*(S02+z*(S03+z*(S04+z*S05))))
	r *= x
	z = 0.5*x + r/s
	if sign {
		return -z
	}
	return z
}

// Y1返回第二类的一阶贝塞尔函数。特殊情况是：
// Y1（+Inf）=0 
// Y1（0）=-Inf 
// Y1（x<0）=NaN 
// Y1（NaN）=NaN 
func Y1(x float64) float64 {
	const (
		TwoM54 = 1.0 / (1 << 54)             // 2**-54 0x3c90000000000000 
		Two129 = 1 << 129                    // 2**129 0x480000000000000 
		U00    = -1.96057090646238940668e-01 // 0xBFC91866143CBC8A 
		U01    = 5.04438716639811282616e-02  // 0x3FA9D3C776292CD1 
		U02    = -1.91256895875763547298e-03 // 0xBF5F55E54844F50F 
		U03    = 2.35252600561610495928e-05  // 0x3EF8AB038FA6B88E 
		V04    = 1.66559246207992079114e-11  // 0x3DB25039DACA772A 
	)
	// 特殊情况
	switch {
	case x < 0 || IsNaN(x):
		return NaN()
	case IsInf(x, 1):
		return 0
	case x == 0:
		return Inf(-1)
	}

	if x >= 2 {
		s, c := Sincos(x)
		ss := -s - c
		cc := s - c

		// 确保x+x不会溢出
		if x < MaxFloat64/2 {
			z := Cos(x + x)
			if s*c > 0 {
				cc = z / ss
			} else {
				ss = z / cc
			}
		}
		// y1（x）=sqrt（2/（pi*x））*（p1（x）*sin x0+q1（x）*cos（x0））
		// sin（x0）=sin（x）cos（3pi/4）-cos（x）sin（3pi/4）
		// /=1/sqrt（2）*（cos（x）+sin x）
		// 为了避免取消，请使用

		var z float64
		if x > Two129 {
			z = (1 / SqrtPi) * ss / Sqrt(x)
		} else {
			u := pone(x)
			v := qone(x)
			z = (1 / SqrtPi) * (u*ss + v*cc) / Sqrt(x)
		}
		return z
	}
	if x <= TwoM54 { // x<2**-54 
		return -(2 / Pi) / x
	}
	z := x * x
	u := U00 + z*(U01+z*(U02+z*(U03+z*U04)))
	v := 1 + z*(V00+z*(V01+z*(V02+z*(V03+z*V04))))
	return x*(u/v) + (2/Pi)*(J1(x)*Log(x)-1/x)
}

// x>=8时，pone的渐近展开式为
// 1+15/128 s**2-4725/2**15 s**4-。。。，式中s=1/x.
// 我们用
// pone（x）=1+（R/s）
// 式中R=pr0+pr1*s**2+pr2*s**4+pr5*s**10 
// s=1+ps0*s**2+…+4.6）

var p1R8 = [6]float64{
}
var p1S8 = [5]float64{
}

var p1R5 = [6]float64{
}
var p1S5 = [5]float64{
}

var p1R3 = [6]float64{
	1.17186865567253592491e-01, // 0x3FBDFFF55B21D17B 
	3.93297750033315640650e+00, // 0x400F76BCE85EAD8A 
	3.51194035591636932736e+01, // 0x40418F489DA6D129 
	9.10550110750781271918e+01, // 0x4056C3854D2C1837 
}
var p1S3 = [5]float64{
	3.47913095001251519989e+01, // 0x40416541349C 
}

var p1R2 = [6]float64{
	1.07710830106873743082e-07, // ，
}
var p1S2 = [5]float64{
	2.32276469057162813669e+02, // 0x406D08D8D5A2DBD9 
	1.17679373287147100768e+02, // 0x405D6B7ADA1884A9 
	8.36463893371618283368e+00, // 0x4020BAB1F44E5192 
}

func pone(x float64) float64 {
	var p *[6]float64
	var q *[5]float64
	if x >= 8 {
		p = &p1R8
		q = &p1S8
	} else if x >= 4.5454 {
		p = &p1R5
		q = &p1S5
	} else if x >= 2.8571 {
		p = &p1R3
		q = &p1S3
	} else if x >= 2 {
		p = &p1R2
		q = &p1S2
	}
	z := 1 / (x * x)
	r := p[0] + z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))))
	s := 1.0 + z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))))
	return 1 + r/s
}

// 对于x>=8，qone的渐近展开式是
// 我们用
// qone（x）=s*（0.375+（R/s））
// 其中R=qr1*s**2+qr2*s**4+qr5*s**10 
// s=1+qs1*s**2+…+6*6*s*6*6*s*6*12 

var q1R8 = [6]float64{
}
var q1S8 = [6]float64{
}

var q1R5 = [6]float64{
	-8.05644828123936029840e+00, // 0xC0201CE6CA03AD4B 
	-1.83669607474888380239e+02, // 0xC066F56D6CA7B9B0 
	-1.37319376065508163265e+03, // 0xC09574C66931734F 
	-2.61244440453215656817e+03, // 0xC0A468E388FDA79D 
}
var q1S5 = [6]float64{
	1.99179873460485964642e+03,  // 0x40F1F31E77BF839 
	1.74684851924908907677e+04,  // 0x40D10F1F0D64CE29 
	4.98514270910352279316e+04,  // 0x40E8576DAABAD197 
	2.79480751638918118260e+04,  // 0x40DB4B04B04CF7C364B 
	-4.71918354795128470869e+03, // 0xC0B26F2EFCFFA004 
}

var q1R3 = [6]float64{
}
var q1S3 = [6]float64{
	-1.35201191444307340817e+02, // 0 0 0 0 0 606060E6060606060606 70290A310 10 10 F ABCDEFF的10 ABCDBBBBBBBBBBBB6 6 6级BCDCDCDCDCDFFFFFFFFFFBBBCDCDCDCDFFFFFFBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBCDG
}

var q1R2 = [6]float64{
}
var q1S2 = [6]float64{
	2.95333629060523854548e+01,  // 0x403D888A78AE64FF 
	2.52981549982190529136e+02,  // 0x406F9F68DB821CBA 
	7.57502834868645436472e+02,  // 0x4087AC05CE49A0F7 
	-4.95949898822628210127e+00, // 0xC013D686E71BE86B 
}

func qone(x float64) float64 {
	var p, q *[6]float64
	if x >= 8 {
		p = &q1R8
		q = &q1S8
	} else if x >= 4.5454 {
		p = &q1R5
		q = &q1S5
	} else if x >= 2.8571 {
		p = &q1R3
		q = &q1S3
	} else if x >= 2 {
		p = &q1R2
		q = &q1S2
	}
	z := 1 / (x * x)
	r := p[0] + z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))))
	s := 1 + z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))))
	return (0.375 + r/s) / x
}
